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I.  Introduction 

In  this  report,  multichannel,  one-dimensional  linear  prediction  is 
considered  where  the  zero-mean  observation  processes  are  assumed  to  be 
generated  by  a  multivariate  autoregressive  model  with  a  white  noise  driving 
vector  u(n)  that  is  partially  correlated  over  the  channels.  It  is  proposed  here  that 
this  white  noise  driving  vector  can  be  modeled  in  terms  of  a  temporally  and 

'spatially'  (ie.,  over  the  channels)  causal  vector  z(n)  uncorrelated  across  channels 
such  that  u(n)=Ly  z(n)  where  is  unit  diagonal  lower  triangular.  A  one¬ 
dimensional  multichannel  prediction  error  filter  of  order  p  is  initially  considered 

H  H 

with  coefficients  of  the  form  Ap(k),  k=0,l,...,p  where  Ap  (k)  is  the  kth  matrix 


coefficient  obtained  from  the  multichannel  normal  equation.  The  output  error 
vector  £p(n)  of  this  filter  has  a  non-diagonal  forward  error  covariance  matrix 

H 


PR££(0)=[Zf]p  which  is  decomposed  by  an  Lyp^Y  Wp  decomposition.  It  is  shown 


-1  -1  -1  H 

that  Ly^  is  an  estimate  of  and  that  a  filter  with  coefficients  LypAp(k), 


k=0,l,...,p  yields  an  output  yp(n)  which  equals  z(n)  when  the  coefficient  estimates 

are  minimum  mean-square  error  (MMSE).  Thus,  a  two  stage  filtering  operatio. 

is  implemented  first  to  obtain  the  coefficients  Ap(k),  k=0,l,...,p  and  [Xf]p 

H  .  -1  . 

followed  by  the  Ly^Dy^Lyp  decomposition  of  [Sf]p  to  obtain  Ly^.  Finally,  it  is 

-1  H 

shown  that  the  coefficients  of  the  form  LypAp(k),  k=0,l,...,p  and  p=0,l,...,P  are 


-1  H 

contained  in  the  block  rows  of  where  Lj^Dj^Lj^  is  a  decomposition  of  the 


multichannel  correlation  matrix  R^x  written  in  index  ordered  form,  L|^  is  unit 

diagonal  lower  triangular  and  is  strictly  diagonal  (ie.,  with  1x1  diagonal 

-1 

elements).  Thus,  the  pth  block  row  of  L|^  contains  the  coefficients  of  the  pth 
order  multichannel  prediction  error  filter.  This  result  implies  that  a  set  of 


1 


-1 


H. 


coefficients  equivalent  to  those  of  Ly  Ap(k),  k=0,l,...,p 


could  be  obtained  in  a 


single  stage  of  computation. 

In  [1,14],  a  model-based  multichannel  likelihood  ratio  is  presented  which 
considers  the  detection  of  a  signal  vector  in  an  additive  non-white  noise  vector 
disturbance  process.  The  implementation  of  this  detection  scheme  utilizes 
prediction  error  filters  to  whiten  the  observation  processes  in  time  and  across 
channels  and  provides  a  practical  framework  to  utilize  the  coefficient  estimates. 
In  addition,  the  whitening  operation  in  time  and  across  channels  provides  a 
generalization  of  the  innovations-based  likelihood  ratio  suggested  in  [1 1,12]  to  the 
mutichannel  case.  As  noted  in  [13],  this  approach  provides  a  generalization  of 
several  detection  problems  as  well  as  a  practical  implementation  architecture. 

It  is  anticipated  that  recognition  of  the  result  presented  here  will  motivate 
the  development  of  computationally  efficient  methods  to  obtain  the  multichannel 
prediction  error  coefficients.  Furthermore,  the  similarity  of  the  processing 
methods  discussed  here  to  those  used  in  two  dimensional  image  processing  is 
noted  in  [9]. 
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II.  Autoregressive  Process  Models 

In  this  analysis,  the  observation  processes  are  assumed  to  be  multichannel 
autoregressive.  TTie  order  M  autoregressive  (AR)  model  for  a  Jxl  vector  process 
x(n)  is 

M  „ 

x(n)  =  -  X  AM(k)x(n-k)  +  u(n)  (2.1) 

k:rj 

H 

where  Aivi(k)  is  the  kth  JxJ  matrix  coefficient  for  the  AR  process  and  H  denotes 

the  Hermitian  operation  (ie.,  the  complex  conjugate  transpose  operation).  We 
H 

note  that  A^Ck)  is  expressed  in  terms  of  the  Hermitian  operation  for  notational 

convenience,  but  is  not  necessarily  a  Hermitian  matrix.  The  vector  u(n)  is  a  Jxl 
white  noise  driving  vector  which,  in  general,  has  an  arbitrary  correlation  across 
the  J  channels  so  that 


r  [0]  1^0 

E(u(n)uH(n  - 1)]  =  I  =  o. 

where  Ruu(0)=[2^f]M  ^^e  JxJ  positive  semi-definite  covariance  matrix  of  the 

forward  white  noise  driving  vector  u(n)  and  may  have  off-diagonal  components. 
Since  u(n)  is  uncorrelated  in  time,  but  retains  an  arbitrary  correlation  across 
channels,  then  with  wide-sense  joint  stationarity  of  the  channel  processes  assumed, 
we  can  consider 

u(n)  =  CyV(n)  (2.3) 

where  the  J  x  J  matrix  Cy  is  a  constant  matrix,  and  v(n)  is  a  white  noise  vector 
uncorrelated  both  in  time  and  across  channels  such  that 


E[y(n)v*^(n  - 1)]  = 


I  ^  0 

I  =  0. 


(2.4) 
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The  elements  of  the  diagonal  matrix  Dy  are  the  variance  terms  associated 
with  the  white  noise  driving  term  on  each  channel.  And  so,  from  eq  (2.3)  we  can 
obtain  the  zero-lag  correlation  matrix 

Ruu(O)  =  E  [u(n)uH(n)|  (2.5a) 

=  E  lC„v(n)  YH(n)Cv  1  (2.5b) 

=  C^DyCv .  (2-5c) 

We  could  assume  unit  variance  on  all  elements  of  Dy  without  loss  of 

generality;  ie.,  Dy=I,  the  JxJ  identity  matrix,  so  that  eq(2.5c)  is  the  Cholesky 

H 

decomposition  Ruu(0)=CyCv  when  Cy  is  lower  triangular.  Since  Ruu(O) 

expressed  in  eq(2.5)  is  Hermitiant  ,  positive  semi-definite,  we  could  also  consider 
an  LDLH  decomposition*  of  Ruu(O)  such  that 

Ruu(O)  =  L^DuL”  (2.6) 

where  Lu  is  unit  diagonal  lower  triangular  and  Du  is  strictly  diagonal  with  non¬ 
negative  entries  only.  Assuming  Dy=I,  the  matrices  Ly  and  Cy  are  related  as 

(2.7) 

SO  that 

C,  =  L„D;".  (2.8) 

Solving  eq(2.6)  for  D^,  we  obtain 


^  It  i.s  noted  that  in  general  the  correlation  matrix  Ruu  (0  's  not  Hcrmitian  for  1^0. 

*  The  motivation  for  considering  the  LDL^  decomposition  stems  from  its  relation  to 
linear  prediction  [3|  and  the  desire  to  determine  both  temporal  and  cross-channel 
whitening  with  a  single  stage  of  decomposition. 


4 


Du  =  Lu'Ruu(0)(Lu')”  (2-9a) 

=  E  lL„'u(n)uH(n)(Lu  )”]  (2.9b) 

=  E  [z(n)zH(n)]  (2.9c) 

where 

z(n)  =  Lu^u(n)  (2.10) 

is  a  J  X  1  vector  containing  uncorrelated  elements.  Note  that  from  eqs(2.3),  (2.8) 
and  (2.10),  we  can  obtain 

z(n)  =  Lu  y(n)  (Zlla) 

=  Di'y(n)  (2.11b) 


when  Dy=l.  Thus  z(n)  represents  a  temporally  and  "spatially-causal"  white  noise 
driving  term  to  the  multichannel  AR  model.  Since  is  also  lower  triangular 
unit  diagonal,  it  is  invertible  so  that  from  eq(2.10) 

u(n)  =  LyZ(n).  (2.12) 

Eq  (2.12)  indicates  that  u(n),  originally  defined  in  eq(2.3),  could  identically  be 
generated  by  the  z(n)  process  through  the  transformation  matrix  Ly;  i.e.,  eq(2.1) 

can  be  written  in  the  equivalent  form 
M  „ 

x(n)  =  -  2  A  M(l^)2L(n-k)  +  z(n)  (2.13) 

k=l 

In  chapter  III,  a  two  stage  multichannel  prediction  error  filter  is 
considered  which  uses  estimates  of  the  A^fk)  coefficients  to  obtain  u(n)  in  the 

first  stage  of  filtering  and  an  estimate  of  to  obtain  the  temporally  and 
spatially  uncorrelated  process  z(n)  in  the  second  stage. 
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Finally,  we  note  that  the  relationship  between  the  matrix  coefficients 
H 

Aivi(k),  the  covariance  matrix  [5^flM  of  the  forward  AR  driving  noise  vector  and 
the  correlation  matrix  R^x  expressed  by  the  multichannel  Yule-Walker 
equation;  ie.. 


Am(R„I  =  |0|...|0|| 

(2.14) 

where 

a"  =|IAm(I)Am(2)  ... 

(2.15a) 

H 

l^xx  -  h[in  n,|vi  4n,n-M  1 

(2.15b) 

in.n-M  =  [xT(n)  xT(n-i)...xT(n-M)l 

(2.15c) 

and 

x'r(k)  =  Ixi(k)  X2(k)...xj(k)l. 

(2.15d) 

The  matrix  [Rxxl  l^e  reversed  order  correlation  matrix  of  [Rxxl;  i-e., 
the  correlation  matrix  obtained  with  the  time  order  of  the  vector  in  n-M  from 

eq(2.15c)  reversed  [10]. 
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III.  Multichannel  Prediction  Error  Filtering 

The  output  of  a  multichannel  linear  prediction  error  filter  of  order  p  ic 

expressed  as 

8p(n)  =  x(n)  -  ^(nIn-1)  (3-la) 

=  x.(n)  +  ^  Ap  (k)  x(n-k)  (3.1b) 

k=l 

=  Ap(k)  x(n-k)  (3-lc) 

k=0 

where  Ap  (k)  is  the  kth  matrix  coefficient  of  the  linear  predictor  with  k=l,2,...,p 
and  Ap  (0)=I,  the  JxJ  identity  matrix.  The  subscript  p  distinguishes  the  matrix 

coefficients  and  the  error  output  as  those  from  a  filter  of  order  p.  In  Appendix 
A,  it  is  shown  that  when  the  matrix  coefficients  in  eq(3.1)  satisfy  the  multichannel 
normal  equations, 

E  [£p(n)  £p  H(n-k)]  =  [0]  k>0,  (3-2) 

where  the  output  vector  process  £p(n)  is  the  MMSE  prediction  error.  The 
multichannel  normal  equations  for  a  prediction  error  filter  of  order  p  are 
expressed  as 

Ap  IRm' =  I'” 

where 

a"  =  II  a"(1)  a"(2)  a"(p)) 

H 

[J^xxJ  “  ^l=n,n-p  =n,n-pl  ~  ^^=n-p,n  =n-p,n 

and 


(3.3a) 

(3.3b) 

]  (3.3c) 
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(3.3d) 


|5:flp  =  E|£p(n)£p»(n)|  =  ‘’ReeCO). 

[R^xl  is  the  reversed  order  correlation  matrix  defined  in  the  previous  section. 

We  also  note  that  eqs  (3.3a)  and  (2.14)  are  identical  in  form.  This  equality 

implies  that  the  MMSE  estimate  of  the  AR  observation  process  is  obtained  when 

the  prediction  error  filter  coefficients  are  identically  equal  to  the  AR  process 

coefficients.  Assuming  the  applicability  of  the  AR  process  as  a  representative 

H 

model  of  the  observation  process,  then  as  the  coefficients  Ap  (k)  k=l,2,...,p 
approach  the  AR  coefficients  of  the  model;  ie.,  for  p>M  and 


tl  ,  ,  |AM(k)  k<M 
Ap  (k)  — 

1(0]  k>M 

we  have  from  eqs(2.1)  and  (3.1a) 


g:p(n)  =  u  (n). 


(3.4) 


(3.5) 


f  igure  3.1  shows  the  synthesis  and  analysis  procedure.  Under  the  condition 
that  the  prediction  error  coefficients  are  identical  to  the  coefficients  of  the  AR 
model  process  and  under  the  assumption  that  the  AR  proc'^ss  is  the  exact  model  of 
the  observed  process,  the  prediction  error  filter  outpi  .  ep(n)  is  a  white  noise 

vector  equivalent  to  the  AR  model  white  noise  driving  vector.  However,  it  must 
be  emphasized  that  the  use  of  an  AR  process  with  a  white  noise  driving  function 
is  usually  an  approximate  representation;  i.e.  it  is  not  used  to  describe  the 
underlying  physical  mechanisms  which  give  rise  to  the  random  processes. 
Rather,  it  is  an  approximate  model  representation  for  these  processes.  We  must 

therefore  make  a  distinction  between  the  model  of  the  processes  (synthesis)  and 
the  estimation  process  (analysis)  [10].  In  general,  the  output  £p(n)  of  the  linear 

predictor  is  not  a  white  noise  vector  output  due  to  the  approximate  representation 
of  physical  processes  by  an  AR  model  as  well  as  the  fact  that  we  often  do  not  have 
'a  priori'  knowledge  regarding  the  values  of  the  coefficients  of  this  approximate 
model.  As  a  result,  we  must  estimate  the  filter  coefficients  from  the  observation 
data  as  we  obtain  it.  With  a  limited  amount  of  data,  the  filter  coefficients  are 
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only  estimates  of  the  AR  process  coefficients.  Furthermore,  assumptions  such  as 
the  ergodic  approximation  are  often  made  in  the  estimation  of  process 
parameters.  In  [2],  the  validity  of  this  approximation  is  considered  in  terms  of 
basic  parameters  of  the  observation  processes.  These  results  provide  a  clearer 
understanding  as  well  as  a  means  of  determining  the  processing  requirements 
needed  to  achieve  parameter  estimates  with  high  confidence  limits. 


SYNTHESIS  ANALYSIS 


Figure  3.1 

For  stationary  processes,  these  coefficients  could  be  determined  through 
estimates  of  the  multichannel  correlation  matrix  and  the  Levinson-Wiggins- 
Robinson  algorithm  [4].  Other  methods  proposed  by  Strand-Nuttall  and  Morf- 
Vieira  have  been  developed  [5,6,7,81  with  improved  performance  with  limited 
data.  For  non-stationary  processes,  adaptive  schemes  must  be  considered.  We 
will  address  this  topic  in  a  subsequent  report. 
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At  this  point,  we  note  that  eq(3.5)  resulted  from  the  analysis  procedure  [via 
a  linear  prediction  error  filter  with  coefficients  given  by  eq(3.4)]  of  the  process 

synthesized  by  eq  (2.1).  If  u(n)  is  assumed  to  be  uncorrelated  across  channels, 
the  resulting  ep(n)  is  also  uncorrelated  in  time  and  space  (i.e.,  channels).  In 

general,  however,  u(n)  may  possess  arbitrary  correlation  between  the  J  channel 
elements.  Therefore,  the  vector  EpCn)  will  retain  a  residual  correlation  over  the 

channels  due  to  the  spatial  correlation  of  u(n). 

Since  the  matrix  *’ReE(0)  =  E  [ep(n)ep  (n)]  is  Hermitiant  ,  and  positive 

semi-definite,  we  can  perform  an  LDLH  decomposition*  of  PRee(O)  such  that 

PRee  (0)  =  Lyp  Dy^  L^p.  (3.6) 

Solving  for  Dy 


V  Lyp  Ree(0)(^p)H 

(3.7a) 

=  l'yp  E  [£p(n)  (4p)H 

(3.7b) 

=E  £p(n)  £p(n)  (Lyp)Hl 

(3.7c) 

=E  |yp(n)  yp  (n)l 

(3.7d) 

where 

ip(n)  =  Lyp  £p(n) 

(3.8) 

so  that  the  vector  yp(n)  contains  uncorrelated  elements.  Also, 


^  II  is  noted  that  in  general  the  correlation  matrix  Rpe  (I)  is  not  Hermitian  for 
*  Other  decompositions  could  be  used  such  as  Cholesky  or  unitary  forms.  However,  as 

noted  previously,  the  motivation  for  using  the  LDL^  decomposition  is  to  note  its 
relation  to  linear  prediction. 


n  n' 


(3.9) 


E  l£p(n)£p(n')!  =  [0] 

from  the  orthogonality  principle.  Then,  using  eq(3.8)  to  solve  for  ep(n)  and 


substituting  this  result  in  (3.9),  we  obtain 

E  :^p(n)  :^p  (n')  Lyp]  =  [0] 

n  n' 

(3.10) 

SO  that 

LypE  [  :ip(n)  :^p  (n')  1  Lyp  =  [0] 

n  ^  n'. 

(3.11) 

Finally,  since  Ly^  is  non-singular. 

E  tepCn)  lip  (n’)  ]  =  [0]  n  #  n'. 

(3.12) 

Eq  (3.12)  implies  that  :yp(n)  retains  its  temporal  decorrelation  while  eq(3.7d) 
denotes  its  spatial  decorrelation.  When 


Ap(k) 


|AM(k)  k<M 
l[0]  k>M 


(3.13) 


the  output  of  the  first  filter  stage  converges  toward  u(n),  so  that  eq  (3.8)  becomes 
(noting  that  Re£(0)  =  Ruu((^)  and  the  uniqueness  of  the  LDL^  decomposition) 


-1  -1 

yp(n)  =  £p(n)  =  Ly  u(n) 

-1 

~  ^u  ^u— 
=  z(n) 


(3.14a) 

(3.14b) 

(3.14c) 


where  eq(2.12)  was  used  to  obtain  eq(3.14b).  Thus,  as  the  filter  coefficients 
converge  to  the  coefficients  of  the  AR  process,  :yp(n)  approximates  the  spatially 


1  1 


and  temporally  whitened  process  z(n).  In  addition,  yp(n)  has  been  obtained 

through  a  causal  and  causally  invertible  transformation  of  the  original 

observation  process  x(n);  ie.,  the  input  vector  x(n)  could  be  recovered  through  an 
inverse  filter  operation  on  :yp(n). 

The  estimate  could  be  obtained  by  first  estimating  the  correlation 

matrix  of  PRee(O);  ie.,  P^^eeCO).  An  LU  decomposition  of  this  matrix  would 

-1  .  ^-1 
provide  Ly  .  The  inverse  matrix  Ly  would  then  be  the  required  estimate  . 


Using  eq(3.1b)  in  (3.8),  we  have 


:yp(n)  =  Ly  ep(n)  =  Ly  [  x(n)  +  Ap  (k)  x(n-k)] 


(3.15a) 


Ly  Ap(k)  x(n  -  k). 


(3.15b) 


The  block  diagram  for  the  filter  implementation  of  eq(3.15a)  is  shown  in 

Fig  3.2  as  a  two  stage  process  with  the  first  stage  consisting  of  the  conventional 
prediction  error  filter  used  to  obtain  £p(n)  while  the  second  stage  consists  of  the 

-I 

operation  on  £p(n)  with  Ly^  to  yield  the  temporally  and  spatially  uncorrelated 

process  yp(n).  The  equivalent  filter  realization  of  eq(3.15b)  is  shown  in  Fig  3.2. 

From  this  viewpoint,  the  MMSE  multichannel  filter  output  which  is  both 

-1  H 

temporally  whitened  and  'spatially'  whitened  consists  of  coefficients  Ly^  Ap(k) 

k=0,l,...,p.  In  the  next  section,  it  is  shown  that  these  coefficients  are  contained  in 

-1  H  . 

the  matrix  elements  of  where  is  used  in  the  decomposition  of 

Thus,  the  motivation  for  considering  the  LDlJ^  decomposition  is  based  on 

the  anticipation  of  utilizing  a  single  stage  recursive  procedure  to  obtain  the  filter 
-1  H 

coefficients  1^  Ap  (k)  required  to  estimate  yp(n).  In  [1,14],  the  process  :^p(n)  is 


used  in  the  development  of  a  multichannel  likelihood  ratio  for  the  detection  of 
partially  correlated  random  signals  in  noise. 


(l+d-u)  xLilJ  i-  I  (7,-u)  X  1-  (i-u)  X  I  1  I  (u)  X 


HANNEL  PEF  STRUCTURE 


IV.  Relation  of  Linear  Prediction  to  Matrix  Decomposition 

In  this  section,  it  is  shown  that  the  coefficients  of  linear  prediction  for  a 
multichannel  random  process  are  related  to  the  multichannel  correlation  matrix 
through  triangular  decomposition.  The  discussion  presented  here  represents  a 

generalization  of  the  scalar  case  noted  in  [3]. 

The  JNxJN  covariance  matrix  R^x  for  a  set  of  N  zero-mean,  stationary 

random  Jxl  vector  processes  x(n)  n=l,2,...,N  can  be  expressed  using  the  form  of 
eqs(2.15b)  through  (2.15d)  with  replacing  eq(2.15c)  and  J  indicating  the 

number  of  channels.  Since  Rxx  is  Hermitian  and  positive  semi-definite,  we  can 
express  it  in  terms  of  a  unique  factorization  such  that 

i^xx  ~  i ) 


where  L^  is  a  JNxJN  unit  diagonal  lower  triangular  matrix  and  is  a  real, 

strictly  diagonal  JNxJN  matrix  (ie.,  with  1x1  diagonal  elements).  Solving  eq(4.1) 
for  D0  ,  we  obtain 

Dg  (4.2a) 


=  E[4  S,,Ni”N(4)“]  (4.2b) 

=  E  1  (4.2c) 

where 


contains  uncorrelated  elements  in  time  and  across  channels, 
expanded  form 


fid) 

~x(l) 

fi(2) 

B2I  ^22 

x(2) 

§1.N  - 

fi(3) 

=r 

B31  B32  B33 

x(3) 

-fi(N)- 

_Bni  Bn2  Bn3  ...  Bj^n_ 

-x(N)_ 

Writing  eq(4.3)  in 


(4.4a) 


where 


k  =  1,2,.. .,N 


(4.4b) 


a'^(k)  =  [Pi(k)  P2(k)-  Pj(k)l 

-1 

and  n  is  a  JxJ  matrix  with  1  <  m,n  <  N  and  m>  n.  is  represented  by  the 

matrix  in  eq(4.4a).  The  matrices  „  are  full,  in  general,  when  m>n,  but  are 
lower  triangular,  unit  diagonal  when  m=n  (ie.,  when  B^p  p  lies  along  the  block 

diagonal  of  the  matrix  in  eq(4.4a)].  It  will  be  shown  that  each  Jx(p+l)J  block 

-1 

row  of  the  matrix  Lp  contains  the  prediction  error  filter  coefficients  for  a 

multichannel  linear  predictor  of  order  p  (for  orders  p=0  through  N-1, 
respectively)  which  decorrelate  the  process  in  time  and  across  channels. 

It  is  now  shown  that  these  coefficients  are  equivalent  to  those  obtained  in 
the  two  stage  whitening  process  developed  in  the  previous  section  resulting  in  the 
temporally  and  spatially  decorrelated  error  vector  process,  :y(n).  This  will  be 
achieved  by  considering  a  two  stage  LDL  decomposition  of  ^^r  stationary 

processes,  we  can  write  eq(4. 1 )  in  block  form  as 


Rxx(O)  B.j^j^(-1)  ...  Rx^x('N+1) 
Rxx(l)  Rxx(O)  .•.Rxx(-N+2) 


Rxx(N-l)  Rxx(N-2)  ...  Rxx(O) 


where  the  JxJ  correlation  matrices  on  the  RHS  of  eq(4.5)  are  expressed  as 


(4.5) 


Rxx(l)  =  E[a(n)xH(n  - 1)| 


(4.6) 


and  the  superscript  B  denotes  block  form.  We  note  that  although  each  Rxx(0 

B 

eq  (4.5)  is  not  Hermitian,  Rj^^  is  block  Hermitian  and  positive  semi-definite  so 
that  a  block  LDL  decomposition  of  eq(4.5)  can  be  made.  We  therefore  consider 

B  H 

~  ^ 

where  L^  is  a  lower  triangular  matrix  with  JxJ  identity  matrices  along  the  block 
diagonal  and  is  block  diagonal.  Solving  for  Dg  , 


D£  =  Ls  rL  (Le)” 


(4.8a) 


=  Lg  E  [  (L^f 

=  E  [  Lp  N  (Lp  )  ] 


(4.8b) 


(4.8c) 


~  E  [£i,n£i,n1 


where 


ii.N  “  Eg  ^i,N- 


(4.8d) 


(4.9) 


In  Appendix  B,  the  block  rows  of  are  shown  to  contain  the  multichannel 
coefficients  for  forward  linear  prediction  for  filters  of  order  p=N-l;  ie., 


-1 


Aid)  I 
Ajd)  A2(1) 
A?(3)  a“(2) 


A^d) 


(4.10) 


aIJCp)  a"(p-1)  Ap(p-2)  ...  Ap(I)  1 


From  eq(4.8d) 


Dg  =  E 


'£(1)1  [£»(1)£»(2)...£H(N)] 
£(2) 

£(3) 

.£(N)_ 


%A0) 


where 


%J0) 


(4.11) 


n  =  p  =  N  -1 


H. 


PR^jCO)  =  E  [£p(n)£;,‘(n)l  s  llflgp  n=l,2 . N  p  =  N -1  (4.12) 

is  the  JxJ  Hermitian  covariance  matrix  of  the  error  vector  for  the  forward 
prediction  error  filter  of  order  p. 

We  now  perform  an  LDL^  decomposition  of  Dg  as, 


H 


(4.13) 


where  is  a  JNxJN  unit  diagonal  lower  triangular  matrix  and  Dp  is  strictly 
diagonal  (ie.,  with  1x1  diagonal  elements).  Solving  for  Dp  , 


Dp  =Lj  D£(L^')H.  (4.14) 

From  eq(4.8d) 

Dp  =  L^'  E  (£,.Nfi“Nl  (4' )“ 

=  E[4'£,,N£“n(4)H]  (4.15b) 

where 

3(i,n  ~  £i,n  “  =i,n  (4.16) 


and  the  RHS  of  eq(4.16)  follows  from  (4.9).  We  will  now  expand  the  matrix 
form  of  the  RriS  of  eq(4.I4).  First,  consider  the  form  of  [the  form  of  is 

derived  in  Appendix  B  and  noted  in  eq(4.10)].  Using  the  JxJ  matrices  „  m^ 

-1 

for  the  off-diagonal  block  matrices  of  , 


Gi,l  [0]  [0] 


[0] 


[0] 


Dp  - 


^2,1  G22  [0] 


Gn,1  Gjsj2 


®Rpp(0) 


'Ree(O) 


[0] 

[0] 

Gn,n  - 

Gpi  G2, 

1  •••  Gn,1  1 

[0]  G2,2  •••  Gn^2 

[01  ... 

10)  g"n_ 

'’Ree(0)_J 


(4.17) 


where  p=N-l.  In  order  to  yield  a  strictly  diagonal  matrix  D£  ,  we  must  have 


and 


^m,n  “  all  m  n 


(4.18a) 

(4.18b) 


where  Dy^  is  the  diagonal  matrix  for  the  pth  block  row  of  Dp.  Since  [Gp^.£p+i] 
is  lower  triangular,  unit  diagonal,  and  since  eq(4.18b)  implies  [Gp+i  p+i]£p(n)  is 
spatially  uncorrelated,  then 


-1 

[Gp+pp+il  =  Ly^ 


(4.19) 


-1 

where  Ly  is  lower  triangular,  unit  diagonal.  Eqs(4.18a)  and  (4.18b)  imply  that 
the  first  matrix  in  eq(4.17)  can  be  written  as 


1  9 


(4.20) 


(4.21a) 


(4.21b) 


and  from  eq(4.16) 


I 

At'd)  I  (0| 

A2(2)  Aid)  I 
A?(3)  A3(2)  A?d)  I 
aE!(P)A?(P-I)A?(P-2)  ...  ,\E!d) 


SO  that 


(4.22a) 
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LYjAi(1) 

H  H 

1\,N  =  LyA2(2)  LyA2(1; 


Ly^pCP)  Ap(P-l)  Ly\p(P-2)  ...  Ap(l)  Ly 


-x(N)- 


A 

Lp  Xj 


(4.22b) 

(4.22c) 


~  1 

where  Lp  is  the  matrix  described  in  eq(4.22b).  From  the  equality  of  eqs(4.16) 


and  (4.22c),  it  follows  that 


-1  _  -]  -1 
Lj.  -  Ly  Lg  . 


(4.23) 


-1  -1 
Since  is  unit  diagonal  lower  triangular,  each  JxJ  matrix  Ly^  is  also  unit 

-I  . 

diagonal  lower  triangular  [see  eq(4.20)].  And  so,  Lp  is  unit  diagonal  lower 
triangular.  Substituting  eq(4.8a)  into  (4.14),  we  obtain 


-1  -IB  -1  H  -1  H 

Dc  =  4  4  Rxx  <4  >  '4  > 

- 1  - 1  ■  ^  H  ■  ^  H 

=  Ly  Lg  Rxx  (Lg  )  (Ly  ) 


(4.24a) 


(4.24b) 


Using  eq(4.23)  in  (4.24b), 


(4.25) 


From  the  uniqueness  of  the  LDL*^  decomposition,  it  follows  from 
eqs(4.2a)  and  (4.24)  that 

Lp  -  (4.26a) 

and 

Dp  =  D^.  (4.26b) 

And  so,  Lp  contains  the  multichannel  coefficients  of  linear  prediction  for  filters 

of  order  p=0  to  N-1.  Furthermore,  these  coefficients  can  be  estimated  either 
from  an  LDL^  decomposition  of  ^  stage  operation  as  denoted  by 

eq(4.16). 
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V.  Summary 

In  this  report,  we  discussed  the  relationship  of  multichannel  linear 
prediction  to  triangular  matrix  decomposition.  Initially,  a  two  stage  prediction 
error  filter  is  described  which  whitens  the  observation  processes  both  in  time  and 
across  channels.  This  filtering  operation  is  such  that  it  retains  a  causal  form 
associated  with  linear  prediction  in  both  dimensions.  For  a  filter  of  order  p, 

■  t  H  H 

these  coefficients  are  LYAp(k)  k=0,l,...,p  where  Ap(0)=I,  the  JxJ  identity 

Ip  w 

matrix.  Next,  these  coefficients  are  shown  to  be  equivalent  to  those  obtained 

TI 

through  an  LDL  decomposition  of  the  multichannel  correlation  matrix. 

Specifically,  the  matrix  coefficients  contained  in  the  matrix  of  eq(4.22b)  are 

-1  - 1 

identical  to  those  in  eq(4.4a).  The  block  rows  of  Lp  and  contain  the 

multichannel  coefficients  for  linear  prediction  filters  of  orders  zero  through 
p=N-l  required  to  'whiten'  the  observation  process  both  in  time  and  space  (ie., 
across  channels).  It  was  shown  that  the  prediction  error  process  that  is  both 
temporally  and  spatially  white  can  be  obtained  in  one  stage  using  the  blocks  from 


also  be  obtained  recursively  by  suitably  modifying  standard  fast  algorithms  such 
as  Levinson-Wiggins-Robinson,  Nuttall-Strand  and  Morf-Vieira  algorithms. 

It  is  anticipated  that  recognition  of  the  result  derived  in  this  report  may 
lead  to  alternative  processing  methods  involving  recursive  algorithms  for 
multichannel  signal  processing  applications  as  well  as  for  two  dimensional  image 
processing  [9]. 
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APPENDIX  A 

In  this  section,  we  determine  the  conditions  under  which  ep(n)  as  expressed 

in  eq(3-l)  is  a  temporal  white  noise  MMSE  output  of  a  linear  prediction  error 
filter.  The  linear  prediction  error  of  x(n)  as  defined  in  eq  (3.1a)  is 


£p(n)  =  x(n)  -  i(nln-l) 


(A.l) 


where  x(nln-l)  represents  the  estimated  vector  of  x(n)  using  past  data  values  with 
the  initial  condition  x(llO)  =  0.  Using  a  linear  predictor  with  p  past  values,  we 
define 


In-1)  =  -  Ap(k)x(n-k) 


(A.2) 


where  Ap  (k),  k  =  l,2,...,p  are  JxJ  matrices  representing  the  coefficients  of  the 
linear  predictor.  Substituting  eq  (A-2)  in  (A-1) 


£p(n)  =  x(n)  + 


p (k)x(n-k) 


(A-3a) 


=  i  Ap(k)x. 


P (k)x(n-k) 


(A-3b) 


where  Ap  (0)  =  I.  Let  the  concatenated  column  vector  of  p+1  vectors  (each  of 


dimension  J)  be  defined  as 

-  rvT/r>\  vT/r 


^n,n-p  [2cHn)  x*(n-l)...x‘(n-p). 

Post  multiplying  eq  ( A-3)  by  j^.p  and  taking  the  expected  value 

H  H  H 

E[£(n)in,n-p^  ~  ^f^P:4n,n-p^,n-pl 
H  H 

-  Ap  El^  j^.pl 


(A.4) 


(A. 5a) 


(A.5b) 


-Ap  [Rxxl 


(A.5c) 
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where  Ap  and  the  reversed  order  correlation  matrix  [Rxxl  arc  defined  in 

eq(2.12a)  and  the  discussion  following  (2.12b),  respectively.  Eq  (A. 5)  can  also 
be  written  such  that 

El£p(n)^^  j^.p]  =  E  {£p(n)  [  x^Cn)  x”(n-l)  xH(n-2)  ...  x^(n-p)]  }  (A.6) 

H 

We  now  determine  the  coefficients  of  Ap  subject  to  the  condition  that  £(n) 
is  a  MMSE  vector.  Under  this  condition, 

E  [£p(n)  xH(n-k)]  =  [0]  k>0.  (A.7) 


Eq  (A.7)  is  the  orthogonality  condition  which  states  that  the  error  vector  is 
orthogonal  to  past  observation  values.  Using  this  condition,  eq  (A.6)  becomes 


E|£p(n)in  n.pl  =  {  E  [£p(n)xH(n)]  |0|  10| ...  |0|  } 


From  eq  (A. 3a) 


x(n)  =  £p(n)  -  X  Ap  (k)x(n-k) 
k=l 

so  that  taking  the  Hermitian  transpose 

x^‘(n)  =  £p  (n)  -  ^  x^(n-k)Ap(k). 
k=l 

Using  eq  (A. 10)  in  the  RHS  of  eq  (A. 8) 


(A.8) 


(A.9) 


(A.  10) 


H 


^[%(n)^n,n-p] 


E(£p(n)£p(n)|-  i  E|£p(n)xH(n-k)|Ap(k) 
k=l 


[0][0]...101} 


I0j...l01 


I 


(A.lla) 

(A.llb) 
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where  we  have  again  used  eq  (A. 7)  to  yield  eq  (A.l  1)  and  is  the  forward 

prediction  error  covariance  matrix  for  the  filter  of  order  p.  Combining  eq  (A-5) 
and  (A-1 1),  we  have 

Ap|Rxxl  =  (ISf|p  I0|[0|...|0|).  (A-12) 

Eq  (A.  12)  is  the  multichannel  AR  Yule-Walker  normal  equation  in  augmented 
form.  It  provides  a  set  of  JP  linear  equations  to  solve  for  the  values  of  the  matrix 
coefficients  which  minimize  the  mean  square  error  vector.  Although  eq  (A-12) 
has  often  been  presented  in  the  literature,  the  reversed  order  form  of  the 
correlation  matrix  has  not  often  been  noted  [101.  We  will  utilize  this  feature  in 
Appendix  B. 

We  now  show  that  the  vector  process  £p(n)  is  uncorrelated  in  time.  At  an 
arbitrary  time  (n-l)  where  l>0,  eq  (A- 10)  becomes 

x^(n-l)  =  £p  (n-1)  -  ^  x^(n-k-l)Ap(k).  I  >  0.  (A.l 3) 

k=l 

Using  eq  (A.  13)  in  eq  (A.7) 

E[£p(n)xH(n-l)l  =  Et£p(a)4‘(n-I)l  -  S  E|£p(n)iH(n-k-l)|Ap(k)  =  |0|.  (A.U) 

k=l 

From  eq  (A7),  we  have 

E  [£p(n)xH(n-k-l)  =  (0|  k>0.  (A.15) 

And  so,  the  summation  term  in  eq(A.14)  is  zero  so  that 

E  lEp(n)£p  (n-l)l  =  |0|  l>0.  (A.  16) 

Thus,  the  sequence  of  vector  outputs  from  the  MMSE  prediction  error  filter  are 
orthogonal.  Since  £p(n)  is  a  zero  mean  Gaussian  process,  its  sequence  of  values 

{ep(n))  are  mutually  independent  so  that  {£p(n)}  is  a  white,  Gaussian  noise 

sequence.  However,  these  processes  are,  in  general,  correlated  across  the 
channels. 
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APPENDIX  B 

In  this  appendix,  we  show  that  the  matrix  coefficients  of  multichannel 
linear  prediction  for  a  multichannel  random  process  and  the  prediction  error 
covariance  matrices  are  related  to  the  covariance  matrix  through  a  block 
triangular  decomposition.  The  procedure  is  a  straight  forward  generalization  of 
[3].  In  this  discussion,  we  are  considering  the  first  stage  of  processing  as  noted  in 
sections  III  and  IV.  This  stage  results  in  the  £p(n)  vector  output.  Therefore,  we 

only  consider  the  Ap(k)  matrix  coefficients  here.  A  treatment  which  includes  the 
second  stage  of  processing  using  the  matrix  Ly  to  obtain  ^pCn)  is  developed  in 


Section  IV. 


B 


Recognizing  that  defined  in  eq(4.5)  is  Hermitian  positive  semi-definite. 


we  can  obtain 

R®,  =  l.£  Dp  (4)H  (B.l) 


where  Le  is  lower  block  triangular  with  the  identity  matrix  I  forming  the  block 
diagonal  matrices  and  De  is  a  real  block  diagonal  matrix.  Solving  eq(B.l)  for 
De 


-1  B  ,  -1  „ 

D£=Lj, 

(B.2) 

If  we  consider 

=  1,N  ~  =1,N 

(B.3) 

we  can  easily  show  that 

H 

Dg  —  E-  l§, ]  jsj  £i,N  1 

(B.4) 

where 

eI'n  =  !£'''(  I)£T(2)...eT(N)| 

(B.5a) 

and 

^],N  =  lx.^(l)  x^(2)...x^(N)l 

(B.5b) 
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Since  Lg  is  lower  block  triangular  with  unit  diagonal  elements,  has  the  same 

form  so  that  eq  (B.3)  is  a  causal  and  causally  invertible  transformation  of  the 
data.  We  now  consider  the  normal  equations  for  a  multichannel  predictor  of 
order  p  such  that 

Ap  1R„I  =  [Iflfip  F  (B-6) 

where  [  ^xxlp  is  the  reversed  order  multichannel  covariance  matrix, 

Ap  =|IAp(l)  Ap(2)..,Ap(p)|  (B.7) 


and 

F=  {1 101(01  ...(0)1  (B.8) 

where  I  is  a  JxJ  identity  matrix.  The  vector  of  matrices  Ap  is  the  vector  of 


multichannel  pth  order  linear  prediction  coefficients  and  [Ef]  is  the  Hermitian, 


pth  order,  multichannel  forward  prediction  error  covariance  matrix.  Post 
multiplying  eq  (B.6)  by  Ap,  and  recognizing  that  is  Hermitian,  we  obtain 

fl  Ap\l)  Ap(2)...Ap(p)][R^] 


Ap(I) 

Ap(2) 


L  Ap(p)  J 


=  I^fle. 


(B.9) 


Using  the  relation 
Gb  Gb  =  I 


(B.IO) 


where  Gb  is  the  block  reflection  matrix  (i.e.,  the  square  J(p+1)  x  J(p+1)  matrix 
with  jxJ  identity  matrices  along  the  block  cross  diagonal),  we  have 


II  Ap(l)  Ap(2)...Ap(p)l  Gb  =  lAp(p)...Ap(2)  Ap(I)  I  1  (B.l  la) 


and 
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(B.llb) 


-  f^xx^- 


Using  these  expressions  in  (B.9),  we  have 
H.  H.  . H 


[Ap(p)...Ap  (2)  Ap(l)Il[R^] 


An(p) 


Ap(2) 

Ap(l) 

-  I  - 


(B.12) 


We  now  write  eq  (B.12)  for  p  =  0,  1, 

"I  n " 

a“(1)  I  [0] 

a“(2)A2(1)  I 

A3”(3)A3(2)  A3”(1)  I 


N-1  in  the  combined  form  as 
^xx(^)  •••  B.xx("N+l) 

Rxx(l)  Rxx(O)  ...  Rxx(-N+2) 

Rxx(2)  Rxx(0  Rxx(^) 


Ap(p)Ap(p-l)Ap(p-2)  ...  Ap(l)  I 


Rxx(N-1)Rxx(N-2) 


Rxx(O)  J 


‘  I  A,(1)A2(2)A3(3)...  Ap(p)  " 

■[Slfleo 

I  A2(1)A3(2)  ...  Ap(p-l) 

I2:fle,  [0] 

I  A3(1)  ...  Ap(p-2) 

= 

•••  •••  ••• 

[0] 

-  [0]  I  - 

(B.13) 

Eq  (B.13)  is  of  the  same  form  as  eq  (B.2).  Since  the  causal  decomposition  in 
eq(B.2)  is  unique^  ,  then  Lg  can  be  identified  with  the  lower  triangular  matrix  in 

eq  (B.13).  llius  the  block  rows  of  are  the  multichannel  coefficients  for 


^  The  uniqucnc.ss  of  this  decomposition  is  based  upon  the  specified  block  form 

defined  using  the  concatenated  vector  in  cq(B.5b).  However,  other  block  forms  of 
could  have  been  made  which  still  retain  the  Hcrmitian  property. 
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linear  predictive  fillers  of  orders  zero  through  p=N-l  and  the  block  diagonal 
matrix  elements  of  Dg_  are  the  prediction  error  variances  associated  with  the 

multichannel  filter  orders. 
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